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ABSTRACT 

Radial velocity and transit surveys indicate that solar-type stars bear super-Earths, 
with mass and period up to ~ 20Mq) and a few months, are more common than 
those with Jupiter-mass gas giants. In many cases, these super-Earths are members 
of multiple-planet systems in which their mutual dynamical interaction has influenced 
their formation and evolution. In this paper, we modify an existing numerical popu- 
lation synthesis scheme to take into account protoplanetary embryos' interaction with 
their evolving natal gaseous disk, as well as their close scatterings and resonant inter- 
action with each other. We show that it is possible for a group of compact embryos to 
emerge interior to the ice line, grow, migrate, and congregate into closely-packed con- 
voys which stall in the proximity of their host stars. After the disk-gas depletion, they 
undergo orbit crossing, close scattering, and giant impacts to form multiple rocky Earths 
or super-Earths in non-resonant orbits around ~ 0.1 AU with moderate eccentricities of 
~ 0.01-0.1. We suggest that most refractory super-Earths with period in the range of 
a few days to weeks may have formed through this process. These super-Earths differ 
from Neptune-like ice giants by their compact sizes and lack of a substantial gaseous 
envelope. 

Subject headings: planetary systems: formation - solar system: formation - stars: 
statics 
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Introduction 



In the past decade, more than 400 extrasolar planets have been discovered around nearby 
solar-type stars. Most of these known planets are probably gas giants because they have masses 
(Mp) and densities (pp) comparable to those of Jupiter and Saturn, while their periods (P) ranges 
from a few days to years. When this population is extrapolated to longer periods (more than a 
decades), their ubiquity implies that the fra ction of solar-type sta rs which have gas giant planets 
around them is approximately t]j ~ 10-15% ()Cumming et al.ll2008l ). 



Gas giant planets must formed in gas-rich environment. However, in contrast to the modest 
value of T?.7, the T Tauri progeni tors of solar-type stars are commonly surrounded by protostellar 
disks (jBeckwith &: SargentI 119961) with surfa ce densities comparable to that of the minimum mass 



solar nebula model (MMSN) (lHavashilll98ll ). Based on the conventional paradigm that gas giant 
planets form in frequently-observed, evolving, natal protostellar disks through the condensation of 
grains, coagulation of planetesimals, mergers of protoplanetary embryos, and gas accretion onto 
super-Earth (with masses in excess of a fe w Earth-masses) co res, we const ructed a population- 
synthesis model through a series of papers dlda fc Linl[2004al lbl. boO.i boOSal lbl. , hereafter Papers 
I-V). In these papers, we qualitatively showed that the low fertility rate of gas giants (i.e., the 
modest value of r/j) around solar type stars despite the omnipresence of planetary cradles around 
their progenitor T Tauri stars is due to 1) an inefficient retention of building-block materials in 
protostellar disks and 2) that gas giant planet formation is only possible in relatively-massive 
(super-MMSN) protostellar disks which deplete on time scales (T^ep) longer than ~ 3 Myr. 

The observed mass-period (M p — P) distribution of these planets appear to be non-uniform 
(jWardI Il986l : ICumming et al.l 120081 ). With a set of well-define and carefully-selected sample, po- 
tential observational se lection effects may be quantitatively identified and taken into account 
(jSchlaufman et al.ll2009l ). Analogous to the stellar color-magnitude diagram for galactic and globu- 
lar clusters, genuine oasis and deserts in planets' Mp — P distribution of a controlled samples can be 
used to cast constraints and provide clues to the theory of planet formation and planetary-system 
evolution. 

Preliminary results of our first sets of simulations essentially reproduced known-planets' ob- 
served Mp — P distribution and the fraction of planet-bearing stars as a function of their mass 
and metallicity (Papers II and I II). They have also be en reproduced and confirmed by similar 
population-synthesis a pproaches (IMo rdasini et al.l |2009| ) and hybrid (N-body + ID viscous-disk 
evolution) simulations (jThommes et al..,2008ai ). albeit there remain some differences in the asymp- 
totic Mp — P distribution. Nevertheless the statistical significance of these agreements remain 
unsettled, because 1) existing observational data are acquired with heterogeneous precession and 
diverse observational selection criteria and 2) theoretical frameworks of some critical physical pro- 
cesses remain poorly understood. Ideally, the construction and upgrades of population-synthesis 
models can be used 1) to set quantitative calibrations on important model parameters and 2) to 
extrapolate tests and predictions for subsequent observations. 
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In our previous five papers, quantitative determination of the efficiency of gas giant formation 
around different types of liost stars, many uncertainties remain in the prescriptions adopted for 
these models. They include the magnitude of the following input parameters as defined in paper 
IV: 

1) the surface density distribution of gas (the fiducial magnitude or fg^ and power index qg 
as defined in Eq. [1] below), effective viscosity {v or the dimensionless a parameter) as a function 
of radius r (see Paper V), mass-flow rate (Mdisk) and depletion time scale Tdep of the disk gas (as 
defined Paper I), 

2) the surface density distribution of solid building-block (grains, planetesimals, and protoplanetary 
embryos) material (the fiducial magnitude or fd and power index qd as defined in Eq. [1]), 

3) the scaling coefficients {h^ and hg as defined in Papers II and III) of dust and gas contents in 
the disk associated with a range of the stellar mass (M^) and metallicity ([Fe/H]), 

4) the rate of type-I migration of protoplanetary embryos which is characterized by a parameter 
Ci in comparison with the linear calculation for unperturbed disks (see Paper IV), 

5) the stalling process and location of both type I and II migration, and 

6) the rate of gas accretion onto protoplanetary cores (which is characterized by ki and k2 in Paper 
I)- 

These properties are determined by the not-so-well understood physical processes such as the 
efficiency of turbulent angular momentum transport (items 1 &: 3) and retention of heavy elements 
in the form of grains (items 2 & 3), planetesimals, and embryos (items 4 & 5). With regard to 
item 6, in the sequential accretion hyp othesis, the formation of gas giaii t s must be precede d by the 
emergence of cores with several M^'s (iMizundllQSG : iPollack et al.lll996l : llkoma et al.ll200Cl ). These 
uncertainties are not independent of each other. In particular, the fraction of solar-type stars which 
bear gas giants depends on the inventory of building block material and the retention efficiency of 
cores. In Paper IV, we demonstrated that in disks with power-law Eg and distributions, the 
observed r]j would be achievable only if 1) the disk is more massive than the MMSN model {i.e. 
^ 1)) and 2) type I migration is relatively inefficient (i.e. Ci < 0.03). 

These requirements need not be satisfied throughout the disk provided there are special loca- 
tions in the disk where building-block grains can accumulate and embryos' type I migration may 
be suppressed (see Paper V). In this context, the ice line appears to be a preferred birth place for 
gas giants because it provides a natural barrier against the orbital dec ay of grains (due to hydro 
dynamic drag) and protoplanetary embryos (due to type I migration) (jKretke fc: LinI 120071 ). The 
period distribution of the simulated models which take into account of the preferred birth place 
of gas giants (Paper V) can match well with a conspic uous (2-3 yrs) up-turn i n observed period 
distributions of Jupiter-mass planets provided Ci ~ 0.1 (jSchlaufman et al.ll2009l ). 



Ongoing transit searches and high-precision radial velocity surveys have led to the discovery of 
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many short-period (with P less than a week) gas giants. According to the population synthesis mod- 
els, these planets relocated to their present-day orbit either through type II migration or dynamical 
instability and relaxation. For this scenario, we need to verify: 1) criteria for gas giants' relocation 
shortly after their formation, 2) efficiency of relocation processes, 3) condition for stalling dynami- 
cal evolution, and 4) extent of gas giants' mass loss and orbital decay during the main sequence life 
span of their host stars. Through detailed quantitative comparisons between population-synthesis 
models and observational data, we can calibrate the efficiency of various relevant and competing 
processes. 

In addition to the reproduction of known observational properties (mostly for gas giants), 
population-synthesis models may also be used to make some robust predictions. A particularly 
noticeable feature in the simulated Mp — P distribution is a domain of planetary desert for inter- 
mediate mass (~ 20 — lOOM®) and period (weeks to months) planets with comparable gas and 
solid internal composition. While the actual boundaries of this parameter region depend on the 
values of various input parameters (such as ki, k2, Ci etc), the paucity of planets in this domain is 
robust, mainly due to the combined effects of runaway gas accretion (Paper I) and type I migration 
(Paper IV). This region of parameter space would be filled if the gas accretion onto the cores is 
significantly slowed down by inefficient heat transfer in the infalling envelope or type I migration 
is grossly inhibited everywhere in the case where planetesimal accretion of migrating cores is not 



inhibited (IMordasini et al.ll2009l ). Thus, forthcoming observations can be used to distinguish key 



assumptions associated with these models. 

In the context of quantitative predictions, the simulated planet distribution also highlights 
rich populations of both Earth-mass rocky planets and long-period ice giants. During the early 
epochs of disks' evolution (especially when they undergo repeated FU Ori outbursts), Mdiski '^g 
and Tid are much larger tha n their values in the MMSN model. In principle, embryos can rapidly 



undergo oligarchic growth (jKokubo fc Idal Il998l . |2002| ) and attain isolation masses Mc^iso which 



are sufficiently large (~ 3 — IOMq) to initiate dynamical gas accretion exterior to the ice line 
(Paper I). However, even with a modest amount of type I migration (i.e. Ci ~ 10~^), the embryos 
would migrate towards their host stars before they reach their isolation masses (Paper IV). In inner 
regions at < 0.5AU, the surface density of rocky materials is severely depleted by type I migration 
of embryos that start migration after they acquire mass comparable to their isolation mass (Paper 
IV). 

The fate of these embryos is determined by their retention probability. If they are efficiently 
retained, they would emerge as conspicuous super-Earths. In contrast to gas giants' type II mi- 
gration, embryos' type I migration and retention processes in the proximity of their host stars are 
poorly understood. In Papers I-V, we counted the number of embryos (with different masses) which 
are expected to cross some inhibited inner disk boundary. 



In this paper, we consider various physical processes which may regulate embryos' dynam- 
ical evolution and determine the formation and retention probability of super-Earths at various 
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locations. This follow-up study is timely because recent radial-velocity surveys with high-cadence 
(several time nightly) and precision (down to ~ m/s) have led to the announced discoveries of sev- 
eral super-Earths with masses and periods up to ~ 2OM0 and a few months (Mayor et al. 2009). 
These data also suggest super-Earths may exist around nearly half of all nearby solar-type stars 
and in many cases in the form of non-resonant multiple-planet systems (M. Mayor, S. Udry, private 
communication). This trend, when confirmed with a well-defined controlled sample, would suggest 
there may be many Earth-mass or sub-Earth slowly-growing cores which failed to accrete any sig- 
nificant amount of gas before their natal disk is evaporated. They also highlight the importance of 
dynamical interaction between multiple embryos during their formation and subsequent evolution. 

In order to analyze the asymptotic evolution of multiple-embryo systems, we briefly summarize, 
in §2, some relevant physical processes and introduce a set of generic prescription which can be 
reliably used to quantitative approximate embryos' interaction with each other and their natal 
disks. In §3, we utilize this prescription, along with our previously established algorithms, to carry 
out some case studies for formation of multiple rocky /icy planetary systems around a majority of 
solar type stars which do not bear any gas giant planets. Based on the results of Paper I and IV, we 
suggest these planets formed in relatively low or modest mass disks. These disks generally have gas 
surface density (S^) less than 2-3 times that of the minimum mass solar nebula (MMSN) model. 
We show that prior to gas depletion, dynamically isolated embryos emerge with sub-Earth masses 
and nearly circular orbits in the inner parts of these disks {i.e. interior to the ice line). We refer 
to this first epoch as the embryos-emergence stage. 

During the (second) migration and accumulation stage, these embryos rapidly migrate to 
proximity of their host stars. Under various circumstances, embryos' migration may be stalled in 
the inner regions of disks around classical T Tauri stars. After the gas depletion, the congregated 
embryos perturb each other's orbits. Within the next ~ 100 Myr, dynamical instability induces 
embryos to cross each other's orbits and growth through cohesive collisions. During this (third) 
giant impact stage, super-Earth form with periods between a few days to months and modest 
eccentricities (e ~ 0.01-0.1). As the host stars mature, inelastic and cohesive collisions between 
embryos with modest relative periapse longitudes damp their eccentricities. In the absence of gas 
giants, systems eventually evolve into a (fourth) stabilizing state in which the residual planets' 
crossing timescale becomes larger than the age of the system. We show that these models naturally 
accounts for the origin of a known population of non-resonant multiple Earth/super-Earth systems 
with periods ranging from days to months. The asymptotic mass of some planets may exceed that 
(several Earth masses) needed to initiate efficient gas accretion in a MMSN environment. But 
these planets acquired most of their masses through giant impacts after the gas is depleted from 
their natal disks. In the absence of gas accretion, such super-Earths do not evolve into gas giants. 
Finally, we summarize our results and discuss their implications in §4. 
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2. Embryos interaction with each other and their natal disks 

The main difference between the present investigation and that in papers I-V are: 1) the 
inclusion of resonant capture between embryos during type I migration, and 2) the calculation of 
embryos' orbital and mass evolution after the gas depletion. These effects are important for the 
formation of multiple Earths and super-Earth systems through collisional merger (giant impact) 
events. We describe in this section some modification and additions to our method of treating these 
interactions. 



2.1. Disk models 



Although disk gas accretion rates in the inner regions and the surface density of dust in the 
outermost regions of protostellar disks have been observationally inferred, there are no reliable ob- 
servational constraints on the gas and heavy-element surface density distribution especially near the 
planet forming regions. Theoretical determination of these quantities relies on poorly determined 
magnitude of angular momentum transfer efficiency. 



In light of this uncertainty, we adopt the MMSN model (jHavashi 



1981 



fiducial set of 

initial conditions. In our population-synthesis model, we introduce multiplicative factors {fd and 
fg) to scale the disk surface densities of gas (S^) and planetesimals (S^^) with those of the MMSN 
model. Following Paper IV, we set 




= Sd,ior/ice/d(r/10AU) 
= S3,lo/3(r/10AU)-9^ 

0.32g/cm^ and Sg,io 



(1) 



where normalization factors T,d,io 
S 

(eq. [3]) and 4.2 for r > oice [the latter can be slightly smaller (~ 3.0) (jPollack et al.lll994l )] 



and at lOAU of the MMSN model, and the step function rj] 



75g/cm^ correspond to 1.4 times of 
= 1 inside the i ce lin e at dice 



In Paper V, we considered disk evolution with constant a and the ice line barrier for migration 
(a local maximum is due to a transition in the thickness of MRI active layers across the ice line). 
In this paper, we use a simple power-law disk model in order to focus our attention on the effects 
of dynamical interactions that we now take into consideration in our sequential planet formation 
model. Nevertheless, we specify an inner disk boundary where T,g vanishes and planetesimals' type 
I migration is arrested (see §2.5). 



Neglecting the detailed energy balance in the disk (jChiang &: GoldreichI 119971 : iGaraud &: Lin 



20071^. we ado pt the equilibrium temperature distribution of optically thin disks prescribed by 
Hayashil (jl98ll ) such that 



280 



lAU 



K, 



(2) 



where and Lq are stellar and solar luminosity. In this simple prescription, we set the ice line to 
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be that determined by an equilibrium temperature (eq. [2]) in optically thin disk regions, 

aice = 2.7(L,/L0)1/2au. (3) 



The magnit ude of Oinp may be modified by the local viscous dissipation ()Lecar et al.ll2006l ) and stellar 
irradiation (jChiang &: Goldreichlll997l : lGaraud &: Linll2007l^. The evolution of the ice line may affect 
frequ ency of gas giant planets around high-mass stars (jKennedy et alJ l2006l : [Kennedy fc Kenyon 
20081 ). These potentially important effects will be incorporated in subsequent papers. 

Dependence of disk metallicity is attributed to distribution of fdfl = fgfiW^^'^^^^d, where fdfl 
and fgfi are initial values of fd and fg, respectively. Due to viscous diffusion and photoevaporation, 
fg decreases with time. For simplicity, we adopt 



fg = /9,oexp(-t/rdep), 



(4) 



where Tdep i s disk lifetime (for detailed disc ussion, see Paper IV). The constant a self-similar solution 



obtained bv lLvnden-Bell Pringld ( 19741 ) is expressed by oc r ^ with an asymptotic exponential 



cut-off at radius tq of the maximum viscous couple. In the region at r < tq, decreases uniformly 
independent of r as the exponential decay does, although the time dependence is different. In 
the self-similar solution. Eg at r < tq decays as T.g oc (t/rjep + 1)"^^^. In the exponential decay 
model that we adopt, decays more rapidly as t/r^cp becomes larger at t > rjep- If the effect 
of photoevaporation is taken into account, decays rapidly after it is significantly depleted, so 
that the exponential decay mimics the effect of photoevaporation. We also carried out some runs 
with Eg cx (t/rdep + 1)"^^^, but we did not find any significant difference in the results from 
the exponential decay cases. In the population synthesis models in Papers I-V, a log uniform 
distribution in a range of 10^-10^ yrs was adopted for Tdep. In the present paper, we fix the value 
Tdep to be 3 X 10^ yrs. 



2.2. Prom Oligarchic Growth to Isolation 



On the basis of oligarchic growth model (jKokubo &: Idal Il998l . |2002| ). growth rate of em- 
bryos/cores at any location a and time t in the presence of disk gas, is described by dMc/dt = 
Mc/Tc^acc where, after correcting some typos in Paper IV, 



2.2 X 10^-^-(^^^/^,,V.-V,-^/^ (^) 



Ma 



Me 



yrs, (5) 



Mc is the mass of the embryo (core), q'^ = qa — 3/2 and q'g = Qg — 3/2, and we have adopted the 
mass of the typical field planetesimals to be m = 10^'^g. 

Embryos' gravitational perturbations continually excite the eccentricity of the field planetesi- 
mals. Prior to the gas depletion, planetesimals' eccentricity is also damped by hydrodynamic drag 
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so that embryos' growth is enhanced by their gravitational foc using effect. In this l i mit, the full 



width of feeding zone of an embryo with a mass is given by (jKokubo &: Idalll998l . |2002| ) 



2Mc 



1/3 



Aa, = 10.H = 10^^j a, (6) 

where rn is two-body Hill radius. In situ growth of embryos is quenched when they acquire all 
the heavy-element building block material within their feeding zone. The m aximum mass (called 



isolat ion mass) of immobile embryos in oligarchic growth stage is given by (iKokubo &: Idalll998l . 

2ooi) 

, . 3/2 „3/2 / a \ 3/4-{3gi/2) / \ '"^ 

In Papers I-V, we adopted the above formulae for mass accretion timescale during the oli- 
garchic growth in the presence of disk gas. We also assumed that after the gas depletion embryos' 
and planetesimals' velocity dispersion becomes comparable to surface escape velocity of the most 
massive embryos at that location. This kinematic transition widens embryos' feeding zone which 
enable them to acquire more planetesimals and attain larger asymptotic masses, though at slower 
rates (because their collisional cross section is reduced to their geometric surface area). 

In these previous first-attempt approximations, each embryo was treated separately and their 
mutual interaction was neglected. For example, in Papers IV and V, decreases in due to accretion 
by independent planetary embryos are computed only in the feeding zones centered around their 
original locations. But, in the presence of disk gas, embryos undergo type I migration due to 
imbalance in the tidal torques from outer and inner disks (see ^2.3p . Relocation to regions with 
different surface densities (which is also depleted by previous generation of embryos) would modify 
embryos' growth rate and asymptotic mass. In order to take into account the possibility of growth 
along migration path, we adopted, in our previous papers, a prescription that embryos' growth rate 
is determined by the same value of fd in their original feeding zone until the migrating embryos 
reach the empirically-derived critical radius (adep,mig) within which the residual planetesimals are 
totally emptied by the preceding migrating embryos. This awkward scheme was introduced for 
computational convenience. In the context of gas-giant formation, embryos evolve into dynamically 
isolated cores with nearly circular orbits in gaseous environment. This approximation is adequate 
for the evaluation of their individual growth and migration. It is also reasonable to assume is 
depleted in this inside-out manner without the consideration of any competitive neighbors. 

In the present paper, however, we intend to determine the mass spectrum and kinematic 
distribution of multiple planetary systems around common host stars. Here, we need to consider 
the concurrent disk evolution along with the mass and dynamical evolution of several coexisting 
embryos in evolving disks (see below and ^2.3p . Under some circumstances, it is possible for the 
Erf decline to be a non monotonic function of a. (For example, the emergence of relatively massive 
embryos across the ice line can lead to the formation of "gaps" in the distribution.) In order to 
take this possibility into account, we compute the evolution of distribution due to accretion by 
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all the emerging embryos in a self-consistent manner. 

In our newly modified scheme, the growth and migration of several embryos are integrated 
simultaneously with the evolution of the S^-distribution. We set up linear grids for across 
the disk with typical width of ~ 10"'^ AU. We introduce a population of seed embryos, all with 
an initial mass 10^'^ g (i.e. that of the residual planetesimals) and compute their mass accretion 
rate. In the inner disk region (interior to the ice line), we set the initial separation of these seed 
planetesimals at any given location oq to be the full feeding-zone width (Aoc = lOrn) of embryos 
with local asymptotic isolation mass Mc^iso- We justify this prescription on the basis that in location 
where rc,acc(-^c,iso) < Tdep; any seed embryos with closer initial separations would have merged and 
additional embryos would have formed unless the entire disk is filled with isolated embryos. Thus, 
this choice of seed planetesimals is optimized for computational efficiency, without the loss of 
completeness. 

In the outer disk region where embryos' growth is slow, they are unlikely to attain a local iso- 
lation masses within the life span of their host stars. A more realistic estimate for their asymptotic 
mass is that (Mg) inferred from equation ([5]) for the local after t ~ 1 Gyr. In the same spirit 
of optimizing computational efficiency without introducing incompleteness, we place seed embryos 
with an initial separation which corresponds to Aoc (full feeding zone width) of embryos with mass 
Mg. For intermediate disk region, the initial spacing for the seed embryos is chosen to be Acc for 
the minimum value of Mc^iso and Mg. Provided there is a sufficient supply of seed planetesimals, 
our results do not depend on the choice of their initial spacing. 

We follow the growth of the seed embryos due to planetesimal accretion in accordance with 
Equations ([5|) and The planetesimals' mass accreted by the embryos is uniformly subtracted 
from the grids that are covered by their instantaneous Acc appropriate for their current mass 
and location. (For example, the feeding zone of an Earth-mass planet at lAU extends over 100 
numerical grid points in the 'Sd distribution.) We follow the evolution of 'Sd in each individual grids 
and use their values to evaluate embryos' accretion rate and the strength of dynamical friction from 
the planetesimals. With this treatment, we no longer need to assume some empirical values for 
and adep,mig in the evaluation of migrating embryos' accretion rate. 

Embryos formed through oligarchic growth attain isolation mass and are well-separated (~ 
Attr) from each ot h er. Despite their mutual gravitational perturbation, tidal drag from disk gas 



(lArtymowica 119931 : IWardI Il993l ) and dynamical friction from planetesimals (e.g., IStewart Sz Ida 



2000) are sufficiently strong to generally preserve embryo s' circular orbits. C onsequently, the sta- 



bility of these widely separated embryos is well preserved (llwasaki et al.ll2002l ). In our prescription, 
the process of embryos' oligarchic growth in a gas-rich environment is adequately represented by 
planetesimal accretion onto initially well-separated seed embryos. Nevertheless, embryos' mutual 
perturbation may become important if their orbits evolve and masses grow beyond isolation (see 
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2.3. Type I migration 



The growing embryos exert tidal torq ue on their natal disks. W i th sufficient masses, they 



unde rgo generally-inward type I migration (jGoldreich &: Tremaine 



1982 



Ward 



1986 



Tanaka et al 



2OO2I ). Several numerical simulations have shown that type I migr ation can be retarded in turbulent 



disks (iLaughlin et al 



laminar disks (iKoller et al 



2004 



Nelson 



200J 



2OO5I: iBaruteau &: LinI I2OIOI ) and weakly viscous or inviscid 



^ Li et al. 2009). The effect of corotation torque and horseshoe 



drag can also reduce the rate and cha.nge the direction of migra tion in disks with various and T 
distributions (jMasset &: Casolill2009l : iPaardekooper et al.ll2010l ). But, the efficiency of these effects 
depends on the turbulent induced angular momentum and heat diffusion coefficients. 

A detailed study of the type I migration is beyond the scope of this paper. However, we can 
utilize our population synthesis models to calibrate the relative importance of several competing 
effects. In this paper, we follow the prescription in P aper IV and use the formula for the migration 
rate that was derived through 3D linear calculation (iTanaka et al.ll2002l ): 



1 



'"migl 



a Ci 2.728 + 1.082g„ \aQ. 



5 X 10"^ X lO^^s 



1 



Clfg \M, 



Mr 



VlAuJ 



3/2 



(8) 



yrs. 



where Ci is a free parameter for the retardation of mi gration, reflecting po ssible non-linear effects 
including turbulence (see Paper IV). The expression of lTanaka et al.l (j2002l ) corresponds to Ci = 1, 
and for slower migration, Ci < 1. It also neglects the effect of horseshoe torque associated with 
either Tig and T distribution. The implication of these effects on the formation of planets and 
planetary systems will be addressed in the next paper together with that associated with a migration 
barrier at the inner boundary of global dead zone. 

In this paper, we are primarily interested in planetary systems without gas giants. Equations 
([5]) and ([7]) indicate that on a time scale ~ Tdep, the most massive embryos emerge from relatively 
low-mass disks (with /g = /d ~ 1) attain masses a few times that of the Earth. Although these 
embryos (or equivalently cores) can accrete gas, it forms slowly contracting envelopes with insignifi- 
cant masses (Paper I). For the case studies in §3, we assume that prior to gas depletion, all embryos 
have sub-critical masses ^ 3 — IOMq and neglect their gas accretion. 

Orbital migration relocates embryos to regions with fresh supplies of residual planetesimals. It 
can also lead embryos to regions where planetesimals are severely depleted by the prior accretion 
of other embryos. Equations ([5]) and ([6]) are formulae for oligarchic growth stage. They indicate 
that both Mc and Mc^iso depend on embryos' local S^. In our numerical scheme, we compute the 
instantaneous local Mc and Mc^iso of migrating embryos together with the evolution of (see 
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2.4. Resonant capture 



In intermediate-mass disks (with 1 ^ /g ^ 3), some embryos form with super-Earth masses. 
Even with a ten-fold reduction in efficiency (with Ci ~ 0.1), these massive embryos undergo type-I 
migration over significant distances from the ice fine prior to the gas depletion. Embryos' relocation 
also provides an opportunity to replenish their feeding zone and to grow beyond their initial isolation 
mass. Differential type I migration can reduce (as well as widen) the orbital separation between 
some embryos. 

Embryos with converging orbits may also capture each other into mean- motion resonances. 
After they enter resonances, these embryos have a tender icy to migrate together while maintaining 



the ratio of their semimajor axes. N-body simulations ([McNeil et al. 



2005 



Odhara fc Idal 120091 ') 



illustrate the possible formation of migrating convoys with several resonant embryos. In order 
to construct an utilitarian prescription for resonant capture, we first compute the rate of type 
I migration for each coexisting embryos separately and independently. We identify convergent 
pairs (i,i) from their differential migration speed. Next, we consider dynamical perturbation 
between embryos of closest pairs Neglecting perturbation by other more distant embryos, 

their orbital separation (6 = \ai — aj\) expands impulsively after each conjunction. For nearly 
circular orbits, the expansion of embryos' spacing after an encounter is given by a linear analysis 
(jGoldreich Tremainelll982l : iHasegawa &: Nakazawalll990l ) as 



5b ~ 30 



(9) 



where rn = ((m^ -|- mj)/3M^y/^a. Since encounters occur at every synodic period [T, 
2TTa/{{dQ/da)b) ~ {AiTa/Sb^lyi)], we find the changing rate to be 



syn 



db 



6b 



~ 7 



i-syn 



b 



VK- 



(10) 



In the limit that db/dt becomes comparable to the differential type I migration speed, Av- 



mig 



j, convergent embryos' mutual interaction would be compensated by their relative 
motion (with speeds Vmig,i and fmigj' respectively). In this equilibrium, the separation of two 
convergent embryos would be maintained at 



I \ 1/6 / A \ -1/4 

btr,p ^ 0.16 ' ' ^ ' "^^^ 



Ma 



VK 



(11) 



For computational simplicity, we set Umig,i = fl/''"mig,i where r^ig^i (for an embryos with a mass rrii 
is given by equation ([S]) , 



6trap ^ 4 X 



</4 



rrii + m 



® 



-1/12 



^-1/4^-1/4 



a \99/4/M 



lAU 



M, 



o 



1/4 



(12) 
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In the above expression, dependences of ^trap on various parameters are very weak and its magnitude 
is always ~ 4 — Srn- 

In the construction of a simple analytic approximation, we assume convergent embryos would 
capture each other into their low-order main-motion resonance with a semimajor axis separation 
comparable to fttrap- We note that if the estimated value of 6trap is smaller than 2\/3'^h> resonant 
capture would not actually occur, because embryos' separation is within their fee ding zone and Ab 
(in eg . [9]) is saturated due to the non-linear effect of overlapping resonances(e.g., iNakazawa &: Ida 
19881 ). 



The condition for resonant capture is marginally satisfied for embryos which undergo rapid type 
I migration with Ci ~ 1. In this limit, trapping would be possible if multiple embryos have similar 
masses and migration speed (so that the ir relative migration speed is slow). N-body simulations 
(jMcNeil et al.l l2005l : lOgihara &: Idal |2009| ) show that several embryos do migrate together in such 
cases. However, near the disk inner edge, earliest-generation embryos accumulate as their migration 
is halted there. When subsequent-generation embryos approach these stalled population with their 
full migration speed, orbit crossing may occur because 5trap is re duced below 4 — Srn . After 



these embryos under go dynamical relaxation and cohesive collisions (jTerquem &: Papaloizoul 120071 : 
Ogihara fc Idal l2009l ) , several embryos may eventually survive as merger products in some new 
resonance configuration. 

In our numerical scheme, we assume that in the limit of inefficient type I migration (with Ci ~ 
0.1), resonant trapping always occurs at 5trap = Srn- Even at the disks' inner edge, the separation 
between stalled resonant Earth-mass embryos is comparable to or larger than the total width of 
their feeding zones. Resonant embryos which captured each other along their migration paths 
would spiral towards their host star in lock-step provided their differential tidal torque continues to 
enforce convergent orbital evolution. During subsequent migr ation, resonant embryos' e ccentricities 
are excited due to the conservation of an adiabatic invariance (iMurray fc Dermottll2000l ). But, they 
are also effectively damped by embryos' tidal interaction with the disk gas. Interior to the ice line 
in MMSN-type disks, an equilibrium is established i n which the largest (Earth-mass) reso nant 



embryos retain a small amount of non circular motion (IChambers et al 



199 



Zhou et al.ll2007l l. In 



general, resonant embryos' equilibrium e is less than the ratio of their separation to their semimajor 
axis Ch = b/a and the resonant embryos' orbits remain dynamically separated. 

Based on this consideration, we neglect the possibility of mergers between resonant embryos 
prior to gas depletion regardless whether they are migrating or stalled near the disks' inner edge. 
We assume converging embryos capture each other into their low-order main-motion resonances for 
which their separation is 6trap — Srjj. In our numerical scheme, we monitor the spacing between all 
migrating embryos. When b between any pairs of resonant embryos is reduced below 6trap during 
orbital integration, we compute, for the next time step, their total angular momentum loss to the 
disk due to type I migration. This loss is then redistributed among the resonant embryos such that 
they would migrate together with a fixed spacing b = Srn between them. 



- 13 - 



When resonant embryos migrate to the inner boundary of the disk, they endure a strong 
corotation torque from the disk gas which halts their orbital decay (see below). This migration 
barrier is maintained during the subsequent arrival and resonant capture of additional incoming 
embryos. 



2.5. Halting migration in the stellar proximity. 



In most inner disk regions (within ~ 1 AU), super-Earth embryos undergo rapid, inward, type 
I migration. However, at special locations where the gas surface density or entropy attain local 
maxima, these embryos' type I m igration may be sta l led due to their tidal inter action with the disk 
gas near their co-orbital region (jMasset et al.ll2006l : iPaardekooper et al.l 120101 ). In MMSN disks, 
a barrier against type I migration is located at the inner boundary of the 'dead zone' where gas 
in the disk midplane is totally neutral and not directly affected by magneto-rotational instabilities 
(jKretke &: LinI 120071 ). However, this barrier moves inward and eventually vanish during the late 
stages of disk evolu tion when become sufficiently small for the entire disk to become active 
(jKretke &: Linll2010l ). Near both (outer and inner) boundaries of the dead zone, angular momentum 
transfer by MRI can induce local uniform rotation in the gas flow. The outer p art of these rigidl y 
rotation zone attains super-Keplerian flow which also induces a migration barrier (iKato et al.ll2009l ). 



Stellar magne tic torque can clear out disk gas inside the radius where it is balanced by viscous 



stress in the disk (Konig! 



1991 



). In typical protostellar disks, this interaction induces the central 

P 

m 



stars to co-rotate with the Keplerian frequency at their ii iner edge. Type I m igration may also be 
halted on the outer boundary of a magnetospheric cavity ( Masset et al.ll2006l ). We discuss, in more 



detail below, a powerful halting mechanism at the outer boundary of the cavity. 

In contrast, gas giants induce gap forma tion in their natal disks a nd undergo type II migration 
which is generally inward inside a few AU's (|Lin &: Papaloizoulll985l ). Type II migration is gener- 
ally stalled well inside the magnetospheric cavity. Since gas is severely depleted throughout this 
magnetospheric cavity, tidal torque between disks and planets well inside the cavity vanishes. In 
the limit of weak magnetic field, the size of the cavity (r^) becomes comparable to or smaller than 
the stellar radii and gas gia nts may be halted by their tidal interaction with their rapidly rotating 
host stars (jLin et al.l 119961 ). Although gas giants are shielded and generally not affected by the 
stellar field direct ly, it can potentially induce Roche-lobe overflow to halt their inward migration 
(jLaine et al.ll2008l ). Unless their host stars' magnetic and spin axes are aligned, close-in gas giants 
would be exposed to periodic field modulation and induced current in the upper envelope on their 
night side where the magnetic diffusivity is relatively large. Ohmic dissipation heats the envelope 
to induce their Roche-lobe overflow. In either case, gas giants are likely to park closer to their host 
stars than super-Earths. 



Thermal expansion of rocky embryos is generally negligible. However, rocky embryos gener- 
ally have higher conductivity than the atmosphere of their host stars so that the magnetic flux 
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tube connecting non-embedded planets and their host stars would slip though the envelope of the 
host stars much faster than across the planets. The potential drop across the field lines drives a 
DC current which is analogous to that proposed for the electrodynamics of the lo- Jupiter system 
(jGoldreich &: Lvnden-Belllll969l ). The Ohmic dissipation of this current produces a torque which 
drives the planet's orbit to evolve toward a state of circularization and synchronization with the 
spin of the star. Around slowly (or rapidly) spinning stars, this process can also cause rocky planets 
with periods less than a few days to undergo orbital decay (or expansion/stagnation) within a few 
Myr, which can affe ct the retention effi ciency of short-period super earths. This effect is discussed 
in a separate paper (jLaine k. Linll20ld ). We neglect it in the present paper. 



In Papers I-V, we artificially halted type II migration of gas giants and type I migration of 
rocky/icy embryos at orbital periods of a few days. This prescription was introduced to simulate 
the effects of an inner cavity in protostellar disks around magnetized T Tauri stars. Using this 
prescription, we were able to infer a statistical estimate on the fraction of "hot jupiters" among all 
known gas giant planets. In this paper, we focus our investigation on the stalling mechanism for 
super-Earths' type I migration at the outer radius (r^) of the magnetospheric cavity. Although the 
transitio n radius between th e dead and active regions is also important for longer-period rocky /icy 
planets ([Kretke &: Linll201Cll ). we will consider this more complex effect in the next paper. 



In order to specify a magnitude for rm, we note that the observed s pin periods of young st ars 



show a bimodal distrib ution at peaks at about one week and one day (|Herbst fc Mundtlbonfj ) Q. 



Herbst Mundtl (j2005l ) suggested that for the one-week period stars the stellar magnetic field is 
coupled with the protostellar disk strongly enough to transfer spin angular momentum to the disk 
and open up a cavity, while the one-day period stars do not have a cavity. In order to take into 
account both possibilities, we adopt two different inner boundary conditions: i) with or ii) without 
cavities. 

In models which neglect the presence of magnetospheric cavity, all the embryos that migrate 
to a radius of 3-day orbital period (0.04AU in the case of = IMq) are removed from the sample 
under the assumption they migrate all the way into their host stars without any stoppage. In 
relatively massive disks, early-generation embryos grow rapidly. Type I migration relocates these 
embryos from their birth place and delivers them to their host stars before they can attain isolation 
mass. Nevertheless, subsequent generations of embryos may form out of the residual planetesimals. 

In order to simulate this continuous formation and migration sequence, we construct a pre- 
scription on the formation of next-generation seed embryos. When type I migration of any embryo 
formed at oq has led to the decay of its orbit to O.Soo, we inject a new seed embryo at oq. In 
contrast to the initial mass (10^'^ g) for the first-generation seed embryos, the initial mass assigned 
to this born-again embryo is chosen to be 10~^ times that of its predecessor. If it is larger than the 



^ Recent Spitzer observations also support the bimodal population ( RebuU et al. 20061 : Cieza fc Baliber 2007 ). 
while Corot observation may suggest single peak population, which means that this issue is still controversial. 



- 15 - 



isolation mass estimated from residual planetesimal surface density, the embryo mass is given by 
the isolation mass. (This prescription is introduced to merely represent the mass growth of residual 
planetesimals.) The growth rate and asymptotic mass of this embryo is determined by of the 
local residual planetesimals. For example, its growth would cease if the total mass of planetesimals 
within its feeding zone is severely depleted. Detailed prescriptions on when to introduce any seed 
embryo and its assigned initial mass do not affect the final result. 

In paper IV and V, we show that prior to severe gas depletion, this self-elimination process 
reduces of the residual planetesimals and the asymptotic mass of late-generation embryos. 
During the gas depletion, a population of embryos is retained in the disk. Many of these late- 
generation embryos have sufficiently low masses to avoid extensive type I migration. In the absence 
of any magnetospheric cavity, since all the early-generation embryos are lost to their host stars, the 
mass and spatial distributions of the asymptotically-retained embryos are not sensitive to the initial 
magnitude of S^. In the inner disk regions where growth time scale of embryos at their isolation 
masses is shorter than their migration time scale, they accrete all the planetesimals in their feeding 
zone. After gas depletion, the retained embryos continue to perturb each other's orbits until they 
undergo orbit crossing and collide with each other (see ^2.6p . In the outer disk region, embryos 
are embedded in residual planetesimals in their feeding zones. Through dynamical friction, these 
planetesimals damp embryos' eccentricities even after the disk gas is completely cleared away (?, 
also see §3.2). 

Around strongly magnetized host stars, surface density of the gas T,g vanishes at r < r^, where 
rm is the radius of the outer boundary of the magnetospheric cavity. Outside this cavity. Eg reaches 
a local maximum at r = r^ax- In the zone at Vm < r < Vmax, pressure gradient in the disk tends 
to be positive and angular momentum is transfered from the disk to the isolated embryos through 
their unsaturated corotation resonances. Nev ertheless, embryos may also lost angular momentum 



to the disk through their Lindblad resonances (iTanaka et al.ll2004^ 



migration is quenched when they attain an torque equilibrium (jMasset et al.ll2006l . Paper V). 



Tanaka &: Ward 



20041 '). Embryos 



When a pair of resonant embryos migrate to the i nner disk edge, they experience a much 
stronger positive torque than that on a single embryo (jOgihara et al.l l2010l ). As we indicated 
above, resonant embryos' eccentricities are excited due to the conservation of an adiabatic in- 
variance. Asymmetric eccentricity damping near the inner disk boundary le ads to a net flow of 



2002 



angular momentum fr om the disk to the embryos. In the linear calculations (jTanaka et al. 
Tanaka &: Wardll2004l ) the timescale of the eccentric damping (re) is shorter than that of type I 



migration (Tmia ) by a factor of ~ (cs/uk)^ ~ 0(10 ^). (In this paper, we consider relatively slow 
migration with Ci ~ 0.1 so that Te/rmig is expected to be even smaller.) Thus, angular momentum 
replenishment to the innermost resonant embryo is sufficient to compensate for the loss of it due to 
type I migration torques on all the trapped resonant embryos. T hrough a series of N-bo dy simula- 
tions which include the damping due to disk-embryo interactions, lOgihara &: Idal (|2009l ) confirmed 
that this "eccentricity trap" can indeed stall the migration of a convoy of resonant planets. Based 
on these results, we consider a limiting case in which the innermost embryos are halted beyond the 
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edge of the cavity. 

In the next section, we only present models in which type I migration is stalled at the outer 
boundary of the magnetospheric cavity. Nevertheless, we consider an alternative prescription in 
which we neglect the effect of eccentricity trap at the edge of the magnetospheric cavity. In that 
case, innermost embryos are often forced into the cavity by the torque from outer embryos. Models 
generated with this prescription predict a large population of very short-period planets which do 
not seems to be consistent with the findings of recent radial velocity surveys. 

For computational convenience, we set the edge of the magnetospheric cavity at 0.04 AU. In 
reality, the location of inner edge depends on the stellar magnetic field and gas accretion rate which 
evolve during the disk depletion. We will consider these evolutionary effects in the next paper. 



2.6. Post-oligarchic growth after gas depletion. 

In Papers I-V, we used Equation ([5]) to compute embryos' planetesimal-accretion rate. This 
formula was accurate for the population synthesis of gas giant planets because we need to determine 
the formation and retention of sufficiently massive cores in gas-rich environments (so that these 
cores can accrete gas and evolve into gas giants). In such an environment, cores' eccentricity is 
suppressed and they attain isolation masses rather than engage in giant impacts. 

Although their building-block embryos may have also formed in gas-rich disks, rocky /icy plan- 
ets' final assemblage need not proceed prior to the gas removal. In typical protostellar disks, 
embryos' eccentricity is effectively damped so that their growth is limited by dynamical isolation. 
However, on the time scale of r^cp ~ a few Myr, gas in these disks is depleted by viscous diffusion or 
photoevaporation while residual planetesimals are exhausted by embryos' accretion except in outer 
regions. With a decline in embryos' eccentricity damping efficiency, their orbi ts become dynarnically 



unstable on time scales which increase with both embryos' separation and Eg ( Chambers et al.lll996l : 



Zhou et al.ll2007l ). Embryos' eccentricities increase until their orbits cross (on a crossing time scale 



g) and their growth resumes through giant impacts. At the end of this post-olig archic growth 



the m asses of the largest embryos typically increase by a factor of several to 30 (e.g.. iKokubo et al 



20061 ). Thus, giant impacts essentially determine the asymptotic properties of rocky/icy planets. 



Although we have not considered the dynamical interaction between multiple embryos in Pa- 
pers I-V, an effort was made to approximate the outcome of giant impacts. In our previous pre- 
scription, the possibility of mergers after the gas depletion was simulated with the expansion of 
embryos' feeding zone. In the construction of mass distribution for close-in rocky/icy planets (Pa- 
per V), we also considered two extreme limits, i.e. either all or none of the embryos migrated to 
the stellar proximity merge. 



In order to accomplish the task to simulate the mass, semimajor axis, and eccentricity dis- 
tributions of multiple super-Earths/Earths systems, we need to construct a prescription which 
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approximates the process of embryos' eccentricity excitation and collisions through giant impacts. 
Our previous application of Equation ^ to the determination of mass growth associated with these 
giant impacts are not accurate. Here, we construct, in the Appendix, an improved prescriptions 
for embryos' eccentricity growth in gas free environment and giant impacts. We outline below our 
computational procedures in a sequential order. 

1) We compile a list of both dynamically isolated embryos and orbit-crossing pairs. We evaluate 
the time scale (rcross, see Appendix) for all dynamically isolated embryos' eccentricity to grow until 
they cross the orbits of their closest neighbors. 

2) We identify the pair of non-orbit -crossing embryos with the shortest Td-oss* 

3) After such a time interval has elapsed, we compute the expected statistical changes in their 
eccentricity and semimajor axis. We then identify all other embryos whose orbits this pair would 
cross if these changes were implemented. 

4) For this group of two or more embryos, we implement statistical changes in e and a due to 
repeated close scattering among themselves. 

5) We apply corrections on the magnitude of semimajor axis changes among the participating 
embryos in order to preserve the conservation of total orbital energy. 

6) We identify pairs of impacting embryos based on their statistically weighted collisional proba- 
bility. 

7) Under the assumption that these events lead to cohesion, we adjust both a and e of the merger 
product to satisfy the conservation of orbital angular momentum. 

The search for potentially orbit-cro s sing p airs fstep 2) are also applied to resonant embryos. 
In this context, IXerquem Papaloizoul (j2007l ) and lOgihara Sz Idal ([20091) found through N-body 
simulations that, for rapid migration (with Ci = 1), several embryos remain locked in mutual mean 
motion resonances near the disk inner edge after gas depletion and dynamical relaxation. These 
results indicate that multiple resonant configuration with relatively small number of bodies is stable 
and the formula for Tcross (eq. [E]]) cannot be accurately applied to such configuration. 



However, lOgihara &: Idal (|2009l ) also found that for migration with Ci ^ 0.1, resonant capture 
is more effective and produce convoys of several (up to dozens) resonant embryos over wide regions 
of the disk ranging from its inner edge to radii beyond 0.1 AU. Resonant configuration of these 
resonant convoys is dynamically unstable and they always start orbit crossing and merge with each 
other after disk gas depletion. Since we are concerned with the relatively slow migration (Ci < 1) 
in this paper, we adopt eq. (fTSl) for Tcross even for the resonance-trapped bodies. 

The above procedure is repeatedly applied while the number of residual embryos declines and 
their separation increases with time. Eventually, the magnitude of Tcross for all residual bodies 
exceeds a Gyr age of their host stars. We classify these asymptotic kinematic properties as dynam- 
ical architecture of mature systems. Although these comprehensive procedures are complicated 
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to integrate, each step is based on well-studied celestial mechanics. Other than two empirical 
parameters, there is no need to introduce any arbitrary assumptions. The two parameters are 
also qualitatively inferred from celestial m echanics (see Append ix), albeit their quantitative values 
are calibrated by N-body simulations bv iKokubo et al.l (|2006l ). Thus, this semi-analytic scheme 
minimizes uncertainties in the embryos' dynamical evolution. 



2.7. Comparisons between the semi-analytic and N-body simulations 

In order to assess the validity of our s emi-analytic scheme, we make direct comparisons with the 



results of analogous N-body simulations. iKokubo et al.l (j2006l ) have performed N-body simulations 



for the giant-impact stage from isolated embryos to terrestrial planets. They systematically varied 
the initial disk mass, or equivalently the isolation mass and number of embryos in the radial range 
of 0.5-1.5AU. Since they carried out 20 runs with different initial azimuthal distribution of the 
embryos for each disk mass, the results of these simulations statistically represent asymptotic state 
of these systems. 

With our prescription, we simulate the evolution of embryos systems with similar sets of initial 
masses and radial distributions. Figures [1] show a typical example of embryos' post-oligarchic 
evolution. In this model, we consider systems which initially contain 16 embryos with Mp ~ 0.1- 
O.2M0. These embryos represent oligarchies which have attained isolation masses in a MMSN-type 
disk {i.e., with /^^ o = !)• The thick and thin lines in the lower panel of Figures [1] correspond to 
the semimajor axes and peri/apo-centers of the embryos' orbits. Close scatterings lead to embryos' 
eccentricity excitation and semimajor axis diffusion. Line discontinuities represent merger events 
between nearby embryos. The upper panel shows semimajor axes and eccentricities of planets in 
an asymptotic state. All planets are assumed to have identical internal density and the radii of 
their representative circles are scaled to be proportional to their physical radii. In this model, 
four planets with masses O.llM® , O.44M0, 1.OM0, and 0.76M® survive in stable orbits at 0.46AU, 
0.56AU, 0.85AU, and 1.74AU, respectively 

In order to consider the statistical properties of these embryos systems, we also simulated a 
set of 20 runs for the same initial fd with our numerical prescription, analogous to the previous 
direct N-body simulations. Figures [2] show the averaged mass, semimajor axis, and eccentricity 
of the most massive bodies and their standard deviations in their final state. For coin parisons, 



these figures include the results obtained by N-body simulations of IKokubo et al.l (|2006l ) and our 
semi-analytical model (panels a and b). We also plotted the same quantities for the second most 
massive bodies in panels c and d. 

The results plotted in these figures include sets of initial conditions with fd = 0.3, 1 and 3. 
Embryos' initial masses are given by their isolation masses and initial semimajor axes are distributed 
from 0.55AU to 1.5AU with orbital separations of lOrn. The total mass of the distributed embryos 
is O.72M0, 2.3M®, and 6.66M0 for fd = 0.3,1, and 3, respectively. According to the values of 



fd, number and ma sses of initial embryos are changed. We used the same initial conditions as 
Kokubo et all (|20od l (see their Table 1). 



The figures show very good agreement in planets' asymptotic mass and semimajor axis between 
the N-body simulations and our models for all initial conditions. The dependence of planets' 
asymptotic mass on fa is well reproduced. Although there is a difference in the expectation value 
for e of the second most massive bodies, it is within a standard deviation. Note that our model 
runs several orders of magnitude faster than the direct N-body simulations, so we can incorporate 
the post-oligarchic dynamical interactions into our synthetic planet formation simulations. 



2.8. Transition from marginal metastability to protracted order 

After the embryos' eccentricity is excited to the magnitude which enables them to cross other's 
orbit, their close encounters lead to either elastic scattering or mergers. (For these self- gravitating 
embryos, we neglect the effect of fragmentation and differentiation between volatile and refractory 
materials.) Repeated elastic scattering further increases embryos' eccentricity to a value Cgsc which 
is the ratio of embryos' surface escape velocity and their Keplerian velocity (see Appendix). With 
this eccentricity, the Safronov number is of the order unity and gravitational focusing no longer 
dominates the collisional cross section. 

Orbit crossing events can also lead to direct collisions and merger events. These cohesive 
collisions reduce the number of surviving planets and widen the separation between them. Although 
some merged embryos' orbits may become temporarily isolated, dynamical instability can re-excite 
their eccentricity. If this marginally stable state can be maintained, the number of surviving massive 
embryos would continue to decline while Tcross would become comparable to the age of the system. 

However, as we discuss in Appendix, e of the merged embryos are often significantly smaller 
than those of their progenitors prior to the collision. These merger events generally involve em- 
bryos with limited range of relative longitude of periastrons centered around 180 degrees. For 
Kepleri an motion with e <C 1 , the mass-weighted total Laplace- Runge-Lenz (LRL) vector is con- 



served (iNakazawa &: Idalll988l . ll989l ). Mergers of embryos with anti parallel longitudes of periastron 
generally lead to fractional cancellation in the mass-weighted total LRL and therefore relatively 
small values of e of the merged embryos (for details, see Appendix). Consequently, the asymptotic 
e of the surviving embryos is generally sma ller than eesr.- This e fficient damping is consistent with 



the results of previous N-body simulations (jKokubo et al.l 120061 ). 



With this efficient damping process, many merged embryos' orbits become temporarily isolated. 
Although dynamical instability can re-excite their eccentricity, the crossing time (tctoss) given by 
eq. ()15p often abruptly increases after some merger events. Figures [3] show time evolution of 
the maximum planetary mass (Mmax) (the top panel), total number (the middle panel) and the 
minimum crossing time (the bottom panel) in five independent runs of our semi-analytical models 
with /d = 1 and initial number of the embryos n = 16. After n decreases to 3-4 and Afmax grows 
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to ~ Mq, Tcross abruptly jumps from Tcross ^ 10^ yrs to Tcross ^ 10^^ yrs, so the systems make a 
transition to a dynamically stable state, on Gyr main-sequence life-span of solar-type stars. Since 
these changes in the mass and eccentricity growth rates are comparable to the gas depletion time 
scale in the disk, stochastic merger events can introduce diversity in the extent of type I migration 
of the surviving embryos. 



3. Population Synthesis of Planetary Systems 



With this new scheme, we simulate the formation of rocky and icy planets. The predicted 
mass, size, period, and eccentricity distribut ion of close-in Earths/s uper-Earths that can be used 
to directly compare with observational data (jSchlaufman et al.ll20ld ). 



3.1. Initial conditions 

In Papers I-V, we presented a series of simulated planetary mass-semimajor axis distributions. 
We adopt a range of disk model parameters which represent the observed distribution of disk 
properties and assign them to each model with an appropriate statistical weight. (For example, 
more massive {fg > 1) and long lasting (rjep > 3 x 10^ yrs) disks are adopted less often than 
low-mass disks.) For each disk, although several planets may be generated over time, their mass 
growth and orbital evolution are treated independently. (Evolution of the distribution was 
simulated with a simplified prescription, see ^2.2p . Mergers and giant impacts after gas depletion 
are approximated by planetesimal-accretion formula for disks' outer regions and artificially enforced 
(or neglected) in the stellar proximity. In Papers I-V, embryos' dynamical interaction with each 
other is neglected in our population synthesis models, and planets' asymptotic mass and period 
distributions are obtained from the compilation of many monte carlo simulations. These properties 
only represent systems which contain single planets, regardless of their masses. 

Solar System contains four terrestrial planets, two gas giants and two ice giants. Many ex- 
trasolar planets have known siblings. Perhaps, all planets reside in multiple-planet systems, albeit 
many members may have sufficiently low mass to be below the detection limit. Nevertheless, their 
dynamical interaction may affect the overall kinematic structure of planetary systems. Even in sys- 
tems which contain only one massive gas giant, its dynamical infiuence on other residual embryos 
may affect their asymptotic architecture. 

In this and subsequent papers, we will use our modified prescriptions described in the last sec- 
tion to generate statistically-weighted mass-period distributions for multiple-planet systems. Since 
there are many competing dynamical processes which may affect the outcome of these simulations, 
we adopt a step-by-step approach with increasing degrees of complexity. Here we first carry out 
some case studies. In these preliminary models, we only consider systems around G dwarfs with 
solar mass and metallicity. (Stellar mass and metallicity dependences will be considered in future 
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investigations.) We focus in this paper on the formation of rocky/icy planets in the absence of any 
gas giants and only consider disks with modest masses (with initial 0.3 < /d = < 3). For these 
low-mass systems, planets' weak secular perturbation on each other has been taken into account in 
our treatment of post-oligarchic growth. (In addition to nonlinear close-encounters, gas giants, if 
present, can also induce sweeping secular and resonant perturbations which will also be considered 
in future papers.) 

The initial distributions of and seed planetesimals are specified in ^2.2i We introduce some 
small random fluctuations to the initial locations of the seed planets for different runs. Mass growth 
and orbital evolution of embryos are computed with the prescriptions in ^2.21^2.51 We do not con- 
sider the possibility of giant impacts among embryos in a gas-rich env i ronme nt be cause these events 



must be preceded by embryos' orbit crossing ( §2.6p . Ilwasaki et al.l (120021 ) and iKominami &: Ida 



(|2002l ) showed that orbital crossing can occur only after T,g is sufficiently depleted so that gas can 
no longer suppress embryos' eccentricity excitation. Here, we adopt a necessary condition for orbit 
crossing to be fg/fgfl < 10^^. Since we use fg = /g,o exp(— t/rdep) with Tdep = 3 x 10^ yrs, orbit 
crossing is possible at t > 2 x 10^ yrs. 

Dynamical friction from a planetesimal swarm can also suppress orbit crossing. We calculate 
the total planetesimal masses in the feeding zones of embryos at each timestep. If the total mass in 
their feed zones is larger than any embryos' mass, we suppress their orbit crossing. (These embryos 
are not included in the orbit crossing bodies in step 4 in the procedures described in Appendix 
§2.6p . Since their growth is slow in outer regions, embryos' perturbation on each other is limited 
and they do not undergo orbit crossing in the absence of any gas giants. We assign an eccentricity 
~ rii/a to each embryo which does not undergo orbit crossing (see step 1 in Appendix). 



3.2. Overall evolution with dynamical interactions 

We first consider a disk with a modest initial mass {fdfl = 2) and migration efficiency Ci = 0.1. 
For illustrative purpose, only seed embryos at 0.5-15AU are integrated in this particular model. We 
apply two different sets of inner boundary conditions: 1) embryos' migration is stalled by a cavity 
(see Figures!!]) or 2) they are not stalled at all. (In subsequent papers, we will consider the possibility 
that some embryos may be forced to undergo further inward migration by the perturbation of 
additional embryos which migrated to their outer low-order main-motion resonances.) The bottom 
panels of these two figures show the time evolution of semimajor axes of all embryos and the top 
and middle panels show eccentricities and masses of final planets at t = 10^ yrs. 

In inner regions, embryo growth due to planetesimal accretion and migration are so fast that 
multiple-generation embryos are formed. In the bottom panels, the lines starting at 10^-10^ yrs 
represent 2nd/3rd-generation embryos. Due to their smaller masses, type I migration of 2nd/3rd- 
generation embryos is relatively slow. They are usually captured and shepherded by the mean- 
motion resonance of Ist-generation embryos that have migrated in from outer regions. 
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At t > 2 X 10^ yrs, gas is sufficiently depleted that it can no longer effectively damp embryos' 
eccentricity. Tliereafter, orbit crossing and coagulations between embryos occur on timescales 
of 10^-10^ yrs. In outer regions, embryos grow slowly and they do not migrate over significant 
distances. Due to dynamical friction by a large population of residual planetesimals, these embryos 
do not start orbit crossing even at t > 2 x 10^ yrs. 

In the model with a magnetospheric cavity (Fig. S]), a convoy of migrating embryos (formed 
interior to a ~ 4AU) congregate and park near the disk edge. Their total mass is ~ 20Mq. These 
embryos are captured into the mean-motion resonances of embryos which have arrived in the stellar 
proximity at earlier times. This group of resonant embryos extend from the inner edge of the disk 



to rad ii beyond 0.1 AU. (This phenomenon is observed in the N-body simulation bv lOgihara &: Ida 



(|2009l ) in their slow migration case.) After disk gas depletion, eccentricity excitation by embryos 
interaction with each other is no longer damped by their tidal interaction with the disk. As their 
eccentricity grows, embryos cross each other's orbit. Close encounters break up their mean motion 
resonance and this group of embryos collide with each other to form six planets with mass in the 
range of ~ 1 — lOM^ between 0.03 and 0.5 AU. In the absence of any residual gas, they cannot 
migrate into resonance again. 

Figure [6] shows the time evolution of the scaling factor (fa) for planetesimal surface density due 
to accretion by embryos in the case of the result in Figure [H In Paper IV, we analytically derived 
the radius beyond which this surface density is nearly preserved to be 

yrs. That analytical formula agrees with the numerical result presented here. In outer regions, 
embryos start their migration when their masses are much smaller than their isolation masses, 
so the depletion in there is small. These residual planetesimals provide damping of embryos' 
eccentricities and thereby suppress orbit crossing and giant impacts in their post-oligarchic stage. 
However, in inner regions (< lAU), most of planetesimals have been accreted by 1 Myr. This 
clearing enables the post-oligarchic growth through orbit crossing and coagulations between isolated 
embryos after gas depletion. 

We also consider a model without any migration barrier at the edge of the magnetospheric 
cavity (see Fig. [5]). In this case, many first generation embryos migrate into their host stars. No 
close-in (< 0.1 AU) Earths/super-Earths survives. Two super-Earths with 0.8M® and 7M^ emerge 
at 0.18 AU and 0.4 AU respectively. 



3.3. Planetary growth in the inner and outer disk regions 



For a solar c omposition, the critical core mass to hydrodynami cally sustain gas envelope is 
M:,hydro ~ WMq (jMizunolllQSd : IPollack et al.lll996l : llkoma et al.ll200d) . Its magn itude also depends 
on the planetesimal accretion rate, atmospheric composition (llkoina et al.ll200ffl and the boundary 



condition between protoplanets' atmosphere and their natal disk (llkoma et al 



2001 



In the stellar 

proximity where the disk gas is relatively dense and hot, magnitude of Mc^hydro rnay be smaller 
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(jlkoma et al.l 120011 ). Once the core mass (Mc) exceeds Mc^hydro envelope starts contracting. If 
Mc > Mc^crit ~ severalM®, the contraction timescale may be less than disk lifetime (Paper I), 
although Mc^crit also depends on atmospheric composition. If embryos arrive in the stellar proximity 
with Mp ^ Mffl, they m ay initiate rapid phase of gas accretion and evolve into hot Jupiters 
(jBodenheimer et al.ll200d l. 

However, embryos' growth in inner regions is regulated by their type I migration and char- 
acterized by a two-stage growth (runaway /oligarchic growth and post-oligarchic growth) process. 
Prior to gas depletion, embryos with mass less than 

M_.0..Cr-(4.)"-°(^)'"(^)-(|;)"Me (-) 

grow in situ, i.e., they gain mass faster than they undergo type I migration. (The critical mass for 
resistance against type I migration Mc^max is determined by the condition Tmigi = Src^acc, where a 

1/3 

factor 3 reflects an actual timescale to reach Mc because Tc acc Mc'.) In regions of disks where 
embryos can attain Mc^max before they acquire all the residual planetesimals within their feeding 
zone would migrate to and accumulate near the disk inner edge. Prior to disk gas depletion, 
since these embryos cannot undergo orbit crossing, their growth through cohesive collisions are 
temporarily quenched. 

Despite being surrounded by gas, these stranded embryos cannot evolve into gas giants in 
disks with modest fd,o{< 10) because the magnitude of Mc^max is well below the critical core mass 

crit) for the onset of efficient gas accretion. Close-in embryos do significantly increase their 
masses through giant impacts to magnitude > -Mc^crit after the disk gas depletion, especially in 
disks with fdfl > 2 — 3. But, there would be little residual gas for these embryos to accrete and 
they are likely to evolve into super-Earths rather than the cores of gas giants. Through this process, 
close-in super-Earths may bypass their isolation masses without becoming gas giants. Thus the 
detection of relatively massive compact planets (with Mp > IOM0) does not necessarily imply a 
high magnitude for M^c crit • 

In the outer regions, on the other hand, planetary growth ends in runaway/oligarchic growth 
stage and type I migration is much less effective. Since isolation mass increases with a, embryos 
emerge outside the ice line can acquire Mp > Mc^core and evolve into the cores of gas giants before gas 
depletion (Paper I). As we stated in the introduction, this paper focus on the low-mass disks which 
do not produce gas giants. In fact, in our disk models with = 2, embryos with a few emerge 
at a ~ 3-5AU (see Figs. [Hand [5]). In the next paper, we will consider dynamical perturbation 
induced on the residual embryos by one or more emerging gas giants (we will also take into account 
the effect of an ice-line barrier in more massive disks, see Paper V). We anticipate a large fraction 
of residual long-period planetesimals and embryos may be ejected as they are scattered by one or 
more gas giants. 
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3.4. Assemblage of rocky planets after gas depletion. 

In this section, we consider the detailed evolution of rocky/icy planets after gas depletion. Us- 
ing the prescription presented above, we present planets' asymptotic Mp — a and e — a distributions 
for several sets of simulations with a range of magnitude in Ci and fdfi- Seed embryos are initially 
distributed between 0.2-20AU (see §3.1). Dynamical interaction between multiple embryos is a 
stochastic process which would be inadequately represented by the results of any single set of sim- 
ulations. In order to characterize planets' statistical distribution for each set of model parameters, 
20 totally independent series of random numbers are used to compute their progenitors' dynamical 
interaction. Prom these Monte Carlo simulations, we obtain the mean values of Mp, a, and e as well 
as their standard deviations. In this analysis, computation with our semi-analytical prescription is 
much faster (by several orders of magnitude) than direct N-body simulations. 

The asymptotic (at t = IGyr) Mp — a (left panels) and e — a (right panels) distributions for 
models with fdfl = 3 and Ci = 0, 0.03, and 0.3 are shown in the top, middle, and bottom panels of 
Figures [71 In these models, we consider highly magnetized host stars and impose a magnetic cavity 
in the disk and assume embryos' type I migration is stalled there (see §2.5p . In §3.51 we consider 
the possibility of negligible stellar magnetic field. 

For presentation purpose, we divide the emerging embryos into the close-in (a < O.laice), inner 
(O.loice < a < 0.3aice), outer terrestrial (O.Saice < a < Oice), and icy {aice < a) planet regions. 
We record number of planets in each region, asymptotic (at t = 1 Gyr) mass, semimajor axis, 
and eccentricity of planets in order of mass in each region. The quantities are averaged over the 
most massive planets, the second most massive planets, the third most massive planets, ... in each 
region. In Figures \7\ we plot planets of the averaged number in each region. 

The results in Figures [7] show that typically three terrestrial (rocky) planets emerge in the 
Ci = (no type I migration) model. Since these planets contain all the building blocks interior to 
the ice line, they have masses of a few M^. (In comparison, for the idealized model in Figure El 
lower-mass planets emerge from planetesimals which were initially located within 1.5 AU in a 
MMSN.) These planets attain relatively small orbital eccentricities (~ 0.1). Their corresponding 
velocity is considerably smaller than their surface escape velocity. In §2.2 and Appendix, we suggest 
that merger events generally lead to some degree of energy dissipation. Nevertheless planets' 
asymptotic eccentricities are la rger than current eccen tricities of Venus and Eart h. Dynamical 



fricti on from residual disk gas (jKominami fc Idall2002l ) or residual planetesimals (| O'Brien et al 



20061 ) may further damp the eccentricities. Because we keep track the amount of residual disk gas 



and planetesimals, these effects can be accounted for in subsequent papers. 

According to equation ()13p . Earth-mass rocky embryos can relocate from all radii interior to 
the ice line to the proximity of their host stars prior to the gas depletion even with inefficient type I 
migration. For models with Ci = 0.03, dozens of embryos form interior to the ice line with masses 
~ Mc^max ~ 0-2Mq) and then migrate to the vicinity of their host star. After gas depletion, these 
embryos undergo dynamical relaxation and cohesive collision to assemble into several super-Earths 
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with a < 0.3 AU. Similar results are produced in models with Ci = 0.3 (see middle and lower panels 
of Figures [7]) . Although more efficient type I migration depletes the residual planet-building blocks 
at several AU, the asymptotic planet distribution close to their host stars is not very different 
between models with Ci = 0.03 and Ci = 0.3. In both models, all the migrated embryos are 
halted near the inner edge. Post depletion giant impacts lead to the assemblage of similar mass 
super-Earths, albeit the number of close-in super-Earths appears to increase with Ci. 

Simulation results for models with less massive disks are shown in Figures [8] {fdfi = fg,o = 1) 
and [9] (/d^o = fgfi = 0.3). Equation (fT3|) shows that in disks with similar metallicity {fd,o/fg,o), 
max increases with fg,o. In principle, small embryos can participate in type I migration. But, 
embryos' growth is also limited by dynamical isolation. Equation ^ indicates that the magnitude 

3 /2 

of the isolation mass Mc^iso also increases with /^ q . In disk regions where -A^ciso ^ -^^c,max5 
many embryos emerge without type I migration. In addition. Equation ([8]) indicates that the 
embryos' migration timescale decreases with both their masses and fg. In the limit of inefficient 
(Ci = 0.03) or no (Ci = 0) migration, embryos would not migrate over significant distance if 
their Tmigi > Tdcp. Limited migration reduces the delivery of building block material to the stellar 
proximity. Consequently, the probability of forming short-period super-Earths is an increasing 
function of fdfi- (The emergence of short-period Earth- mass planets from modest to low-mass 
disks requires relatively efficient type I migration with Ci > 0.3). 

In contrast, the retention of embryos near their cradles (just interior to the ice line) promotes 
the formation of rocky planets with intermediate periods (months to years). After gas depletion, 
post-oligarchic growth continues through giant impacts on a time scale which is an increasing 
function of a. In a MMSN-like disk (with f^fi = 1), a system 4-5 planets emerge (on time scale of 
~ 100 Myr) with comparable masses, semimajor axes, and slightly larger eccentricities to those of 
the terrestrial planets in the solar system. Due to the difference in the impact of type I migration, 
the formation probability of potentially habitable planets in this model is actually higher than that 
out of more massive {fdfi = 3) disks (see Figures [7]). 



3.5. Formation of super-Earths around protostars with weak magnetic fields 

In the limit of negligible stellar magnetic field, protostellar disks extend to the stellar surface 
(see §2.5p . In the absence of a magnetospheric cavity, migrating embryos are unable to be halted. 
The importance of this inner boundary condition is highlighted in Figures [10] which represent 
models with the same value of fdfli= 3) as the model in Figures [71 (In this series of simulations, 
the onset of efficient gas accretion is artificially suppressed. In a subsequent paper, we will consider 
the concurrent formation of gas giants and super-Earths.) 

With modest migration efficiencies (Ci = 0.03 and 0.3), embryos formed within ~ 1 — 2 AU 
prior to gas depletion attain mass ~ Mc^max and migrate all the way into their host stars. They 
either accrete residual planetesimals or capture other embryos along their migration paths and clear 
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the inner (< 1 AU) region around their host stars. At a > 2AU, embryos' migration is initiated 
before they reach their isolation masses (see Paper IV). After the gas depletion, multi-generation 
embryos with non-negligible masses are retained. These embryos undergo orbit crossing as their 
eccentricity is excited by their mutual perturbations. Through post-gas-depletion giant impacts, 
they merge into Earth-mass planets. Beyond the ice line, cores at 4-5 AU have masses of ~ lOM^. 
But, these systems do not contain any Earth- mass planets with periods less than a few month. 
This asymptotic dynamical architecture (i.e., Mp — a distribution) is similar to that of the Solar 
System. 



3.6. Production of debris disks 

Observationally, there is no clear correlation between the detection of gas giant planets and 
that of debris disks around common host stars . While the fraction of s olar type stars with known 



gas giants increases with their metallicity (e.g.. iFischer &: Valenti II2005I . also see Paper II), there is 



no analogous correla tion between stellar metallicity and the presence of debris disks around them 



(|Greaves et al.l 120071 ) 



Debris disks are composed of grains with sizes comparable to the infrared wavelengths (a few 
/im to sub- millimeter) . Around mature solar stars, smallest (sub-micron) grains are blown out by 
radiation pressure. Due to Pointing-Robertson drag, modest (sub-mm) gra ins undergo orbital decay 



on ti me scale (a few Myr) shorter than the age of their host stars (e.g., ITakeuchi Artvmowicz 



2001 



). Thus, these grains must be continually generated through collisions of their parent bodies. 



In all the models we have considered here (regardless the magnitude of Ci), embryos growth 
prior to gas depletion is limited by dynamical isolation in the inner disk regions. After gas de- 
pletion, post-oligarchic giant impacts inevitably occur, regardless of the stellar metallicity. Thus, 
the production of "warm" debris dust (with wavelength up to ~ lO/xm) is expected to be common 
around all solar type stars. 

If, in the inner disk region, a large fraction of the total mass in heavy elements is attained by the 
largest embryos, time scale generating a large pool of warm debris particles (rcross) would be longer 
than the grains' clearing time scale but shorter than the age of mature host stars (Gyr). Thus, 
replenishment of detectable grains in debris disks may be a stochastic process and the intensity of 



debris disk signature may vary episodically (jKenvon &: Bromlevll2004l ). 



In outer disk regions, planetesimal growth is incomplete with a large amount of residual plan- 
etesimals (see §3.2). Embryos' eccentricities in this region are likely to be excited by close encounters 
or merger events with other embryos and damped by gas through tidal interaction and residual 
planetesimals through dynamical friction. Unless icy embryos become sufficiently massive to ef- 
ficiently accretion gas (more massive than Mc^crit), their surface escape velocity (~ 10 km s^^) 
is generally smaller than Keplerian velocity at a < 10 AU. Consequently, close-encounters with 
embryos do not lead to the ejection of residual planetesimals. In this limit, planetesimals' velocity 
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dispersion is excited by embryos' repeated gravitational perturbation and damped by inelastic col- 
lisions among themselves. When an equilibrium is established, planetesimals' velocity dispersion 
becomes comparable to the surface escape speed of the planetesimal population which carries most 
of the masses. 



If these planetesimals have sizes larger than a few km's (as in the Kuio er Belt. (IPan &: Sari 



2005 . e.g.,), their velocity dispersion would exceed that for disruptive collisions (jStewart Leinhardt 
2009 ). Debris of such collisions would produce "cold " grains which reprocess stellar radiation and 
generate far infrared excess radiation (jWvattl l2008l ) . Since the collision time scale for numerous 
km-size planetesimals is expected to be shorter than the grain depletion time scale, the signature 
of outer debris disk is expected to be maintained at a steady level. 

In relatively massive disks where one or more embryos can attain critical mass for the onset 
of efficient gas accretion prior to the depletion of the disk gas, the formation of gas giants can 
scatter residual planetesimals to large distances. For example, in the Solar system, while residual 
planetesimals near Uranus and Neptune were scattered by them to the Kuiper Belt region or the 
Oort's cloud, those i n the vicinity of Jupiter a nd Saturn were either ejected or scattered to the 
distant Oort's clouds (jPuncan k. Levison 1119971 ). The clearing of residual planetesimals (as parent 
bodies) would suppress the signature of cold debris disks around stars with gas giants. Thus, 
around host stars with relatively low metallicity, the lack of gas giants does not necessarily imply 
a lower detection probability for debris disks. 



3.7. Formation of close-in super-Earths 

The results in previous sections indicate that a small amount of planet-disk tidal interaction can 
lead to significant type I migration for embryos more massive than Mc^max- In sufficiently massive 
disks (with f^fi > 1), a convoy of embryos stall from outer edge of magnetospheric cavity (specified 
to be 0.04 AU in our models) to locations beyond 0.1 AU. (Around strongly magnetized stars or 
in less massive disks, embryos may be stalled at larger distances from their host stars). These 
embryos merge through giant impacts after gas depletion. Since only a fraction of embryos' energy 
is dissipated, the resultant semi major axis of the merger products is expected to be comparable 
to that of their progenitor embryos. 

Close encounters between embryos generally do not lead to mean motion resonance. In the 
absence of residual disk gas, these merger products do not undergo any further orbital decay and 
generally remain out of mean motion resonance with each other. These expectations are consistent 
with the simulated asymptotic Mp — a distribution (see Figs. [7] and (H). 

The velocity dispersion of the residual embryos is a fraction of their surface escape speed. In 
the stellar proximity, it is much smaller than the local Keplerian speed. In comparison with Earth- 
mass planets at around lAU, the simulated eccentricities of close-in super-Earths are relatively 
small (in the range of 0.01-0.1). In subsequent papers, we will consider the possibility of excitation 
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of much larger eccentricity by secular resonan ces of distant gas giants (e.g.. iMardling &: Linll2004l : 
Nagasawa et al.ll2005l : iThommes et al. Il2008bl ). 



Many emerging planets have masses larger than ~ Mc^crit which is the critical mass for cores 
to evolve into gas giants within a few Myr. Since these planets acquired most of their mass after 
the gas is depleted or in gas-free cavities, they cannot accrete a substantial gaseous envelope. The 
atmosphere of their progenitor embryos may also be ejected during giant impacts. Despite their 
challenges in attaining primary atmospheres, the emerged super-Earths may attain a metal-rich 
secondary atmosphere through outgassing. However, such an atmosphere would contribute to a 
small fraction of these planets' total mass. 



3.8. Frequency of habitable planets 

In disks with Eg and Y,^ comparable to those of MMSN, embryos' isolation mass at the hab- 
itable zone (~ 1 AU) is considerably smaller than that of the Earth. Nevertheless, terrestrial 
planets attain most of their asymptotic masses, after the gas depletion, through giant collisions 
and merger events as their velocity dispersion increases until they bypass their dynamical isolation. 
This conjecture is support ed by the late Earth's formation ep och (~ 30 — 60 Myr as inferred from 



radioactive isotopes (e.g., lYin et al.ll2002l : iKleine et al.ll2002l )) which appears to persist well after 



the disk gas is depleted (on a few Myr time scale). 

In Papers I-IV, this post-oligarchic growth has been taken into account with a simple prescrip- 
tion. In that approach, the width of embryos' feeding zone (~ 2eesca) is assumed to expand with 
their eccentricity and embryos' asymptotic mass is estimated to be (Paper I) 

This approximation naturally reproduces the dependences of the outcome on the magnitude of /^ q 
and a. In comparison the results simulated with our new prescription, the above expression for 
Me^iso slightly underestimates it (by a factor of ~ 2 or so). 

Figure 4 in Paper V shows that a fraction of solar type stars harboring habitable planets with 
Mp ~ 0.3 - lOM® and a ~ 0.75 - 1.8AU, ry® ~ 20% for Ci < 0.03 and ~ 10% for Ci < 0.3. These 
values do not significantly change by using the present new model. Thus, with our new prescription, 
we verify that the magnitude of rj^ is adequately reproduced by the simulations presented in Paper 
V, using Eq. ^ 



4. Summary and Discussions 



In this paper, we introduce a new prescription for our population synthesis models. In this 
revision, we incorporate, for the first time, a semi analytic prescription for close scatterings and giant 
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impacts among solid planetary embryos (for detailed prescriptions, see Appendix). Our analytical 
model quantitatively reproduces statistics of the asymptotic distributions of mass, semimajor axi s 



and orbital eccentricities of bone fide planets obtained by N-body simulations (IKokubo et al.ll2006l ) , 
as shown in §2.7. 

We also take into account the effects of resonant trapping by embryos during their migration 
and when they are stalled by migration barriers (§2.4, 2.5). With these tools, we are able to 
simulate the statistical distributions of solid planets' asymptotic mass, semimajor axis and orbital 
eccentricities, in systems without any giant planets. 

These simulations are particularly relevant to the recently discovered population of close-in 
super-Earths. Our results indicate that the emergence of super-Earths in the proximity of their 
host stars proceeds through three stages: 



• In a gas rich environment, embryos accrete planetesimals in their feeding zones and become 
dynamically isolated (runway/oligarchic growth). 

• Embryos undergo orbital decay and accumulate close to their host stars (migration/stall). 

• After disk gas depletion, embryos perturb and excite eccentricities of each other. Eventually, 
their orbits cross and they grow significantly beyond the isolation mass through collisions 
among them (post-oligarchic giant impacts). 



During the second and third stages, sub-Earth embryos evolve into systems of non-resonant, 
multiple super-Earth planets through the following mechanisms: 



• The coupled effects of type I migration of embryos, termination of the migration near the 
inner disk edge, and resonant trapping between embryos produce a resonantly trapped convoy 
of embryos that extends beyond 0.1 AU for the the cavity inner boundary condition (§3.2). 

• In the presence of the disk gas, because their orbits are stabilized by efficient tidal e-damping, 
embryos' individual masses are similar to the critical mass to resist against type I migration 
that is < They do not go into runaway gas accretion phase. 

• Through close scattering and giant impacts after disk gas depletion, the resonantly clustered 
embryos form multiple Earths or super-Earths in non-resonant orbits around ~ 0.1 AU in the 
disks with masses comparable to or a few times larger than the MMSN {fdfl ~ 1-3). 

• Merger products attain non-resonant orbits and their orbit crossing timescale is longer than 
the typical age of their host stars (~ 10 Gyr). Consequently, their asymptotic dynamical 
architecture is very stable. 

• Although the Earths/super-Earths suffer close scattering, their eccentricities are as low as 
0.01-0.1 due to efficient collision damping and large local Keplerian velocity. 
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• Although the asymptotic mass of some super-Earths is larger than that required to initiate 
efficient gas accretion, they do not evolve into gas giants because their final assemblage occur 
only after the residual disk gas is depleted. 

In disks with fd^ ^0.3, the initial isolation mass of embryos interior to the ice line is consid- 
erably less than that of the Earth. Type I migration of sub-Earth embryos is less efficient and the 
amount of embryos accumulation in the stellar proximity is limited. Consequently, close-in planets, 
if formed at all, tend to have masses less than that of the Earth. 

For a range of disk accretion rate (comparable to that observed in typical T Tauri stars), 
migration o f embryos formed in the outer regions of protostellar disks is stalled by a barrier near 



the ice line (iKretke &: Linll2007l . Paper V). Although the isolation mass at large distances from their 
host stars is substantially larger than that of the Earth, embryos' growth time scale is also very 
long. Embryos may attain masses (mostly made of icy material) comparable to that of Uranus and 
Neptune but they would not be able to completely clear their feeding zone (even on a time scale of 
a few Gyr). Consequently, their eccentricity is effectively damped by dynamical friction from the 
residual planetesimals even after the disk gas is depleted. 

In our simulations, we found a relatively large fraction of stars bear ice giants, albeit many of 
these stars do not contain any gas giants. While ice giants tend not to cross each other's orbits, they, 
nonetheless, scatter residual planetesimals to large distances (beyond 10-100 All's, analogous to the 
Kuiper Belt) and continuously supply parent bodies of debris disks which emit persistent (over Gyr 
time scale) mid and far IR reprocessed radiation (equivalent to cold zodiacal light). In the region 
interior to the ice line, embryos' growth is suppressed by dynamical isolation rather than collisional 
time scales and a few super-Earth oligarchies emerge while residual planetesimals between them are 
mostly acquired by them. In this region, post-oligarchic collisions occur stochastically and produc e 



dusty fragments which provide sources for the NIR reprocessed radiation (jKenvon Bromlevll2004l ) . 
The occurrence of inner debris disks only weakly depends on the host stars' metallicity and mass 
through the asymptotic mass and spacing of the embryos. 

The cases presented here are for disks with modest mass. In relatively massive {fdfl ^2 — 3) 
or metal-rich disks, one or more gas giants may form (Paper I, IV). We have not considered here 
systems which contain gas giants. The effects of their secular and resonant perturbation will be 
considered in the next paper. Here we draw to the attention that our Solar System does not have 
any planetary body inside Mercury's orbit at 0.4AU. Yet it must have formed from a disk with 
fdfl comparable to or larger than unity (i.e., at least in a MMSN). Many extrasolar planetary 
systems also do not display any sign of close-in super-Earths. We suggest this dichotomy is due 
to a dispersion in the sizes of an inner disk cavity. The extent of this disk region depends on the 
diverse disk accretion rate and strength of host stars' magnetic field during their T Tauri phase of 
evolution. 

In §2.5, we suggested that the observed bimodal distribution of spin periods of young stars 
with peaks at about one week and one day respectively may be produced by the different strength 



- 31 - 



of the coupling (jHerbst fc Mundtll2005l ) . If the couphng is stronger than some critical value (around 
stars with a few kilogauss fields), the stellar spin angular momentum would be transfered to the 
disk through the stellar dipole magnetic field and the ionized disk gas would be accreted through 
the channel flow along the magnetic field lines (e.g-. IShu et al.lll994l ). Inside the radius (~ 0.06 — 0.1 
AU) where disk gas and stellar spin co-rotate (with period of a w eek), the disk w ould lose angular 
momentum and be truncated with a magnetospheric cavity (e.g.. lShu et al.lll994l ). 



This magnetospheric cavity may not exist or be confined to much smaller radii around stars 
with relatively weak magnetic fields. Inefficiency of angular momentum transfer would result in 
much faster stellar spin (with periods of a day). In this limit, migrating embryos would be halted 
(if at all) so close to their host stars that the star-planet tidal interaction would lead to further 
orbital decay and disruption. We simulated such a possibility with models in which no migration 
barrier is imposed at the disk inner boundary (see §3.2 and 3.4). We find that regions interior to 
~ 0.3 are eff'ectively cleared by early generation of migrating embryos. 

First-generation embryos may have formed in the primordial solar nebula when the infant Sun 
had a relatively weak field. Extrasolar planetary systems without close-in super-Earths may also 
have formed around weakly magnetized T Taur i stars. The streng th of magnetic coupling may 
also evolve with the accretion rate in the disk (iKretke &: Linll2010l ). In a subsequent paper, we 
will address the magnetic coupling process that would play a key role in the origin of diversity of 
close-in super-Earths or Earths. 
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Appendix. Prescription for eccentricity excitation and merging of embryos in 

post-oligarchic stage 

We describe here prescriptions for embryos' eccentricity excitations and merging process. Be- 
cause giant impacts proceed after the depletion of disk gas and planetesimals, we consider only 
mutual gravitational interactions among the embryos and neglect damping forces from disk gas 
and dynamical friction from planetesimals. 

The order of our procedures is: 

1) set up initial conditions, 

2) evaluate, for all pairs of embryos, the crossing time scale (rcross) over which sufficiently large 
eccentricities are excited (by their mutual secular perturbations) so that their orbits would cross 
each other, 

3) evaluate the amount of eccentricity excitation (e) and semimajor axis change {5a) of embryos 
pairs which are destined to cross each other's orbits next, 

4) determine whether these changes may lead to the participating embryos to undergo secondary 
orbit crossings with any other neighboring embryos, 

5) compute changes in the e and a for embryos which undergo secondary orbit crossings, 

6) make adjustment in a for all embryos participated in orbit crossing so that their total orbital 
energy is conserved, 

7) identify a pair of orbit-crossing embryos which are most likely to physically collide with each 
other, 

8) create a merged embryo, and 

9) repeat step 2, until Tcross of embryo pairs exceed the system age. 

The detailed prescriptions for each step are as follows: 

1. The initial mass (rrij) and semimajor axis (aj) of jth embryo are determined by oligarchic 
growth model (see section 2.6). The eccentricity (cj) of dynamically isolated embryos are 
assigned randomly by a Rayleigh distribution with mean values of their Hill eccentricities 
(defined as (mj/3M*)"'^/^ where is the host star mass). Embryos' initial eccentricities are 
of the order of 0.001-0.01. Since eccentricities are significantly excited by embryos' mutual 
interaction after the onset of orbit crossing, the initial small values of ej do not affect the 
results. 



2. Given mj,aj, and ej, the crossing time (timescales on which orbital crossing sta rts) of each 
pair of the embryos are calculated, following the fitting formula obtained bv lZhou et al. 
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(|2007l ) with some shght modifications 

log 



A + Blog 



2.3ru, 

i{= ^aittj) of the pair, h 



where Tk is Keplerian time at the mean a{ 
rrii )/3M^y/^a, and 

A = -2 + eo-0.271og;U, 

B = 18.7 + 1. Hog /X- (16.8 + 1.2 log /i)eo, 



1 [Ci + ej)a 



(15) 

{{rrii + 

(16) 



{mi + mj)/2 



3. The pair with the shortest orbit crossing time (t*,.^^^) is assumed to undergo close encounters 
before any other pairs. During these close encounters, ei nbryos' excited velocity dispersion is 
limited by surface escape velocity of the perturber Vescj (|Safronovlll969l : lAarseth et al.lll993l ). 
The corresponding eccentricity is 

1/3 / N 1/6 



V, 



-esc,«j 



esc, J J 
VK 



^2G{rai + mj)/{Ri + Rj 



^GAU/a 



~ 0.28 



rrii + f^t 



P 



a \i/2 

A/m / V3ecm-^/ VlAU/ ' 

(17) 

where vk is Keplerian velocity, and Rj is the physical radius calculated from rrij with bulk 
density p of 3 gcm~^ inside the ice line and 1 gcm~^ outside it. The eccentricity c hange is 



partitioned according to the mass of interacting bodies ([Nakazawa k, Ida 



1988 



eccentricity distribution is relaxed to a Rayleigh distribution (jida &: Making 
the individual excited eccentricities are given by 

nii 



1989) and the 



19921 ). Thus, 



rrii + "ij 



(18) 



where 7^ is a random number produced by a Rayleigh distribution with the mean value of 
unity. Changes in the semimajor axis associated with the eccentricity excitation are assumed 
to be zizejUj. Sign of the change is chosen to ensure that orbital separation between the pair 
is widened after their orbit crossing. 

4. Energy and angular momentum exchanges during the close encounters between any pair of 
orbit-crossing embryos may induce them to cross the orbits of other nearby embryos. We 
classify all the affected embryos by their overlapping (relative to the modified pair) radial 
excursion (epicycle amplitude) into closely interacting groups. These groups generally contain 
several (> 3) embryos. We assume that the most massive embryo in each group dominates the 
outcome of successive close encounters within the group. The characteristic post-encounter 
eccentricity (ej) for all but the most-massive embryos is given by equation (jlSp in which i is 
used to represent the most massive embryos in the group. We also assume the post-encounter 
eccentricity of the most-massive embryo is determined mostly by the second-most-massive 
embryo. 
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In general, successive close encounters lead to radial diffusion of the orbit-crossing embryos. 
In our prescription, changes in embryos' semimajor axis are given by values chosen Wjejaj. 
The statistical weight Wj (which has a skewed distribution within a limited range between 
-1 and 1) is introduced to reflects the kinematic distribution of the orbit crossing embryos. 
If the total mass (Min) of all the encountering embryos in the group inside a particular 
embryo's orbit is more than that (Mout) outside it, its probability of orbital expansion would 
be larger than that of orbital decay. In order to simulate this effect, we introduce a parameter 
fj = Cj(Min - Mout)/(Min + Mout) and set Wj = fj + {1 - \ fj\) * Rj where Rj is a random 
number with a uniform distribution in the range of [—1,1]. With the results of N-body 
simulations, we also calibrate the magnitude of a constant Cj ~ 2/3. 

Finally, close-encounters also excite embryos' inclination . Magnitude of the orbital inclina- 



tions in radian may be about half of that of eccentricities (jlda &: Makindll992l ). However, the 
orbital inclinations do not directly affect the condition for orbit crossing, we are not concerned 
with embryos' inclinations in this paper. 

5. After the relative changes in the semimajor axes of all the orbit crossing embryos have been 
evaluated, we renormalize their semimajor axes by a constant numerical factor so that the 
total orbital energy of all the embryos in the group is conserved. In general, the modification 
of the semimajor axes introduced by this renormalization procedure is less than a percent 
order. This simple algorithm does introduce artificial changes in the total angular momentum 
of the encountering embryos. 

Although, with a more complicated non-uniform renormalization factor, we can also conserve 
the total angular momentum of all the embryos in the group, we opted to limit the degrees of 
freedom at the expense of a small change in the total angular momentum of the system. After 
multiple close encounters and merger events, the total departure from angular momentum 
conservation of the entire system is limited to within 10% in typical models. 

6. Close encounters also leads to physical collisions. We assume all physical collisions are cohesive 
and lead to merger events. Among the embryos in the orbit crossing group, a colliding pair 
{k, I) is chosen. However, after the impulsive momentum changes during repeated close- 
encounters not every chosen pairs from this group would retain overlapping orbits. If the 
radial excursions of any chosen pair do not overlap, another pair would be chosen until all 
available candidates are exhausted. (For example, after a close encounter between embryos 1 
and 2, the former may cross the orbits of embryos 3 and 4 whereas the latter may cross the 
orbits of embryos 3 and 5. If coUisional candidates are chosen to be embryos 1 and 5, this 
pair would be rejected.) 

In order to take into account various factors which contributes to the collisional probabilities, 
we select the appropriate embryos with a statistically weighted random number. N-body 
simulations show that collision frequency decreases with embryos' semimajor axes (because 
Keplerian period is longer and embryos' spatial density is lower at larger semimajor axis). 
We adopt a weighted collisional probability which is proportional to a^'^. 
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7. A merged embryo is created from each colliding pair. Since we neglect fragmentation and 
rebound, the merged embryo acquires the total mass of the colliding pair. In the direct N- 
body simulations, the center of inertia and the total momentum are conserved after each 
cohesive collision. In our prescription, we do not follow the exact location and kinematics of 
the merger product. Instead, we use conservation laws to determine orbital elements of the 
merged embryos. 

In principle, a small fraction of orbital angular momentum is transformed into spin angular 
momentum after each collision. However, the spin angular momentum is usually negligible. 
Total kinetic energy of the colliding embryos is also not conserved, because the collision 
dissipates kinetic energy of the relative motion and the merging changes binding energy. 
However, if their velocity dispersions are smaller than their two-body surface escape velocity 
{vesc,ki), the collisional dissipation energy is approximately equal to the change in the binding 
energy. In practice, the total orbital energy of the colliding pair is approximately conserved 
during the collision as well as the total angular momentum. 

Assuming the conservation of orbital energy, the semimajor axis of the merged body (aki) is 
given by 

ruk + mi _m^^ma ^^^^ 
iki «fc ai 

In the case of e ^ 1, it is difficult to accurately evaluate the eccentricity of the merged body 
from the angular momentum conservation. Here, we assume that after each collision, the 
mass-weighted total Laplace-Runge-Lenz vector is conserved. This assumpt ion holds for col- 



lisions betw een embryos with nearly circular Keplerian motion [i.e. e ^ 1) (iNakazawa &: Ida 



19881 . 11983 ). In this limit, the eccentricity of the merged body (cfc/) is given by 

inik + mi)eki cos Wki = nikek cos Wk + miei cos wi, ^^q^ 
(nik + mi)eki sinwki = nikek sinwk + miei sintu;, 

where w's are longitudes of periastron at the collision. 

Due to the secular perturbation between orbit-crossing embryos, and wi generally precess 
independently over all angles. In principle, both longitudes can be chosen with uniformly 
distributed random phase (between and 27r). However, even with overlapping regions of 
radial excursion, the orbits of a pair of eccentric embryos would only cross each other if the 
relative angle between their longitudes of periastron (6 = tUfe — wi) is within some limited 
range. For example, eccentric embryos would not cross each other's orbits if their longitudes 
of periastron are aligned with ~ 0. In contrast, orbit crossing generally occurs between 
embryos with 6 in some limited range centered at 180 degrees, provided e^afc and ejaj are 
larger than | Ofc — ai\. 

In our prescription, we randomly generate a value of ro/ between and 27r. We assume 
secular perturbation generally induces the longitudes of periastron of any pairs of embryos 
with overlapping orbits to precess and circulate so that their 6 can always enter into the range 
(AO) which allows them to collide. We first compute the magnitude of with independent 
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Keplerian orbits for the neighboring embryos. We then take into account that embryos' 
gravitational perturbation on each other tends to shghtly broaden this range. The value of 6 
at the collision is randomly chosen from a value within the range A^. The value of Wk is then 
specified by wi and 6. Since the preferred value for is ~ vr, the summation in the r. h. s. of 
eq. (pO|) usually results in cancellation, so that e^i is often significantly smaller than and 

8. For updating the system time, we need to take into consideration 1) the time interval for 
neighboring embryos to cross each other's orbits (rcross = ''"cross equation [T5]) and that 
required for embryos to collide with each other Tcoiude- Before each collision, embryos with 
overlapping orbits undergo repeated close scattering as 6 circulates in and out the range of 
AO. The magnitude of Tcoiiide for embryos with overlapping orbits is determined by their area 
filling factor (see Paper I). N-body simulations indicate that at a fraction of AU, Tcoiudc ~ 
10^'^ — IO^'^Tk for Earth-size embryos. 

For pairs with moderately large initial separation (6 > Srn), Tcross ^ "^collide and it is adequate 
to update the evolution of the system from t to t + t*^oss- Pairs with small initial separation 
(6 < Sfh) are closely packed with wide range of A9. Since it takes longer for them to collide 
than to cross each other's orbit, we update the evolution with a randomly generated Tcoiude 
which has a uniform distribution in the range (10^'^ — 10^'^)Tk. Thus, the new system time 
is given by t + r*ross for T*ross > IO^'^Tk and t + T*Qiiide otherwise. 

After some collisional events, it is possible for the merger product to become dynamically 
isolated from all other embryos in the group. N-body simulations also indicate that following 
other cohesive collisions, several additional merger events may follow in rapid succession. For 
example, in a group of 4 orbit-crossing embryos, embryosl and 2 may merge first. The newly 
form embryo may attain a new orbit which continues to overlap those of embryos 3 and 4. In 
that case, it is possible for the newly merged embryo to subsequently coagulate with embryos 
3 or 4 or both. In order to take these possibilities into account, we repeat step 2 to 7 for all 
group members after each merger event, until their t*^.^^^ (including the collisional lag time) 
has exceeded the integration time. 



In step 7, we indicate that the merger product attains an eccentricity (cki) which is often 
significantly smaller either than and ei. This apparent eccentricity damping is due to the effect 
that merger events can only occur for a limited range between the longitude of periastron of the 
colliding embryos. Consequently, embryos' asymptotic eccentricities are usu ally smaller than eps r. 



This inference is consistent with the results of previous N-body simulations (jKokubo et al.ll2006l ). 



In the prescription presented here, the only free empirical parameters are i) the weighted 
changes in the semimajor axis in step 4 and ii) the broadened range of AO in step 6. All the 
other procedures are based on fundamental dynamics although some of their efficiencies are ap- 
proximated. The weighted a changes and the broadening of the phase range are also qualitatively 
based on celestial mechanics. The quantitative parameters for these processes are determined by 
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comparison with the results obtained by N-body simulations (jKokubo et al.ll2006l ). Note that the 
same parameters are applied for all models in the present paper and we do not individually carry 
out fine tuning for different disk conditions. 
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Fig. 1. — An example of evolution of embryos in the post-oligarchic stage calculated by the 
semi-analytical model. Initially, 16 embryos are placed in the disk with = 1. Thick and thin 
lines in the lower panel represent the evolution of embryos' semimajor axes and peri/apo-centers 
respectively. Discontinuities in lines represent merger events for embryos with others. The upper 
panel shows semimajor axes and eccentricities of final planets. The radii of circles are proportional 
to physical sizes. 
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Fig. 2. — The averaged quantities of 20 models obtained by our semi-analytical model are 
represented w i th fil led squares. For comparison purpose, the results of N-body simulations of 
Kokubo et al.l ([200a) are represented with open circles. The bars indicate standard deviations, (a) 
orbital eccentricity and mass of the most massive bodies, (b) their semimajor axis and mass, (c) 
orbital eccentricity and mass of the second most massive bodies, and (d) their semimajor axis and 
mass, the second most massive bodies (panels c and d). Results for three sets of disk parameters 
ifd = 0.3, 1 and 3) are plotted. The average masses increase with fa- 
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Fig. 3. — Time evolution of the mass of the most massive embryo (Mmax) is plotted on the top 
panel. The total number of residual embryos is plotted in the middle panel. The minimum orbit 
crossing time for any pairs of embryos is plotted in the bottom panel. Five sets of models (generated 
with different random number seeds) for /rf = 1 disks are presented here. Each model starts with 
16 embryos at the onset of the simulation. 
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Fig. 4. — Time evolution of semimajor axes of all embryos is plotted in the bottom panel. The 
asymptotic (at 1 Gyr) eccentricities and masses of all "final" planets are plotted in the top and 
middle panels. In this model, we set /^^ = 2. We also assume the presence of a disk cavity is 
sufficient to stall the inward migration of all embryos in the proximity of their host stars. 
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Fig. 5. — The same as figure [5] except we did not impose a migration barrier near the host star. 
In our prescription, we adopt the no-cavity condition for the inner boundary. 
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Fig. 6. — The time evolution of planetesimal surface density due to their accretion by embryos. 
Same model parameters are used as for the result in Figure H The distributions at t = 0, 10^, 10^ 
and 10^ yrs are plotted (the lines at t = and 10^ yrs almost overlap with each other). 
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Fig. 7. — Simulated planetary systems for Ci = 0, 0.03, and 0.3 are shown in the top, middle, and 
bottom panels respectively. These results are generated from models with fdfi = 3. A migration 
barrier is imposed in accordance with the the cavity condition. The left and right panels show 
asymptotic masses and eccentricities of bone fide planets at i = 1 Gyr. The mean values of 20 runs 
are expressed by the symbols with bars of standard deviations. 
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Fig. 8. — The same as Figures [7] except for fd = 1. 
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Fig. 9. — The same as Figures [7] except for fd = 0.3. 
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Fig. 10.— 



The same as Figures [7] except for the non-cavity condition. 



